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distribution  of  mud  based  on  echo-soundings.  The  dissipation  of  waves  by  a  non-rigid  bottom  is 
represented  in  the  wave  model  by  treating  the  mud  layer  as  a  viscous  fluid  Applied  for  431  time 
periods,  the  model  without  this  type  of  dissipation  has  a  strong  tendency  to  overpredict  nearshore  wave 
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the  wave  model,  yields  the  observed  waveheights. 
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1.  Introduction 

Wind-generated  surface  waves  in  shallow  and  intermediate 
depths  generate  pressure  variations  at  the  seabed  with  spatial  and 
temporal  scales  corresponding  to  the  wavelength  and  wave 
period.  In  the  case  of  a  non-rigid  bottom,  such  as  mud,  the 
pressure  variations  can  result  in  motion  of  the  water/seabed 
interface.  Work  is  being  done  to  generate  this  motion,  and  thus 
energy  is  lost  from  the  wind-waves.  In  the  case  of  a  muddy 
bottom,  the  motion  in  the  seabed  is  subsequently  damped, 
predominately  by  viscosity. 

Methods  exist  for  estimating  the  damping  of  water  waves  by 
viscous  mud.  An  early  effort  was  made  by  Cade  (1958),  using  an 
assumption  of  shallow  water.  Dalrymple  and  Liu  ( 1978)  developed 
a  more  general  method  without  using  this  assumption;  further, 
their  method  accounts  for  viscosity  in  the  water,  rather  than  just 
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the  mud  layer.  Ng  (2000)  proposed  a  numerical  simplification  of 
the  Dalrymple  and  Liu  (1978)  calculation,  using  an  assumption  of 
a  thin  mud  layer.  Such  treatments  of  non-rigid  seafloor  as  a 
viscous  fluid  do  have  limitations:  mud  can  also  exhibit  viscoelas¬ 
tic  or  plastic  behavior:  see  Hsiao  and  Shemdin  (1980);  Jiang  and 
Mehta  (1995.  1996);  Zhang  and  Ng  (2006);  Mei  and  Liu  (1987), 
and  references  therein. 

Treatment  of  damping  by  wave-bottom  interaction  within  an 
analytical  wave  model  requires  the  a  priori  assumption  that 
unrepresented  processes  (refraction,  shoaling,  wind  effects,  break 
ing,  etc.)  are  small.  For  verification  with  field  data,  this  assump¬ 
tion  means  that  test  cases  must  be  very  carefully  selected,  with 
most  data  sets  being  unsuitable  Treatment  within  a  numerical 
wave  model  greatly  improves  this  situation  since  these  processes 
can  be  efficiently  incorporated.  One  such  model  is  the  SWAN  wave 
model,  introduced  in  the  1990s  (Holthuijsen  et  al„  1993;  Ris  1997; 
Booij  et  al.,  1999)  for  the  purpose  of  predicting  wave  propagation, 
growth,  and  decay  in  coastal  regions,  and  has  since  seen 
considerable  use  by  scientists  and  engineers.  However,  this  model 
includes  parameterizations  for  attenuation  via  interaction  with  a 
rigid  seafloor  only.  Even  in  early  evaluations  of  the  model,  it  was 
noted  that  the  absence  of  mud-related  dissipation  is  a  major 
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Nomenclature 

Tp  peak  wave  period 

( i  wave  angular  frequency 

D(0)  normalized  directional  distribution  of  wave 
energy 

k  wavenumber 

h  water  layer  thickness  (without  motion) 

(Vo  fluidized  mud  layer  thickness  (without  motion) 
<Vo.<  total  mud  layer  thickness  (without  motion) 
pw  water  density 

mud  density 
V  ratio  y  =  pwfo)m 

vw  kinematic  viscosity  of  water  layer 

vm  kinematic  viscosity  of  mud  layer 

Am  Stokes  boundary  layer  thickness,  Am  ■-  \fvm/a 

{  ratio  C  =  Aml Aw  =  (vm/vw)1/2 

(?  normalized  mud  layer  depth  (Reynolds  number). 

d  =  ^m.ol^m 


k,  imaginary  part  of  wavenumber,  equivalent  to  the 

dissipation  rate 

Dm/vv  or  k,  dissipation  rate  from  viscosity  in  the  mud  and  water 
layers,  in  the  case  of  Ng  and  DL  formulae,  Dni/w  =  k, 
Dm  dissipation  rate  from  viscosity  in  the  mud  layer,  in  the 
case  of  WDGL,  Dm  =  k, 

Dw  dissipation  rate  from  viscosity  in  the  water  layer 
ks  real  part  of  the  wave  number  in  shallow  water,  k5  = 

ff/v'gh 

xa  o.soo  distance  over  which  a  wave  will  be  reduced  to  50%  of 
its  original  amplitude  (see  Eq.  (1)) 
xa  osqo  distance  over  which  a  wave  will  be  reduced  to  10%  of 
its  original  amplitude  (see  Eq.  (1)) 

Sds  spectral  dissipation  rate 

Sbor  spectral  dissipation  rate  from  wave-bottom  interac¬ 

tion,  a  subset  of  Sds 

Smud  spectral  dissipation  rate  from  mud.  a  subset  of  Sho( 

E  spectral  energy  density 

N  spectral  action  density  N  =  E/a 


deficiency,  such  that  in  cases  of  non-rigid  seafloor,  one  must  apply 
unrealistic  bottom  friction  parameters  to  get  the  desired  dissipa¬ 
tion  (e.g.  Dingemans,  1998).  Representation  of  damping  by  mud 
was  recently  introduced  in  SWAN  by  Winterwerp  et  al. 
(2007),  based  on  an  extension  by  De  Wit  (1995)  of  the  Cade 
(1958)  formulation,  though  at  time  of  writing,  it  has  not  been 
incorporated  into  publicly  released  versions  of  the  code  (Hoi  thuijsen 
et  al.,  2006).  Winterwerp  et  al.  (2007)  utilize  laboratory  experi¬ 
ments  with  measured  rheology  by  De  Wit  (1995).  No  prior 
numerical  wave  modeling  study  involves  the  application  of  field 
measurements  of  rheology;  in  Winterwerp  et  al.  (2007),  the  field 
rheology  is  assumed. 

The  primary  objective  of  the  present  study  is  to  use  numerical 
wave  models  to  simulate  wave  dissipation  by  fluid  mud,  utilizing 
field  measurements  of  both  mud  and  wave  conditions.  The  field 
experiment  was  held  at  Cassino  Beach,  Brazil  during  May  June 
2005.  Large  numbers  of  simulations  with  the  SWAN  wave 
model  over  a  35-day  period  are  performed  to  identify  trends 
and  sensitivity  to  physics  of  wave  dissipation  by  viscous  mud  in 
this  application.  Secondary  objectives  are  as  follows:  (i)  to 
compare  two  methods  for  representing  the  dissipation  of  wind¬ 
generated  surface  waves  by  a  viscous  mud  layer:  the  method  of 
Winterwerp  et  al  (2007)  and  that  of  Ng  (2000),  the  latter 
implemented  in  an  experimental  version  of  the  SWAN  model  in 
this  study;  (li)  to  determine  whether  including  dissipation  by  mud 
is  necessary  for  accurately  reproducing  observed  wave  heights; 
(iii)  to,  given  an  estimate  of  rheology  and  mud  distribution 
derived  from  field  measurements,  evaluate  the  skill  of  these 
models  in  predicting  observed  wave  heights;  and  (iv)  to  develop 
and  apply  an  inversion  process  to  determine  the  rheology  and 
mud  distribution  for  which  the  wave  model  will  reproduce  the 
observed  wave  heights. 

Section  2  of  this  manuscript  introduces  the  modeling  platform, 
.SWAN,  as  well  as  the  two  methods  of  representing  dissipation 
by  mud  in  this  platform.  The  methods  are  also  verified  in  this 
section,  including  a  comparison  with  results  using  the  approach 
of  Dalrymple  and  Liu  (1978).  In  Section  3,  the  Cassino  Beach  case 
study  is  introduced.  In  Section  4,  the  two-dimensional  model 
design  is  described,  and  results  are  presented.  In  Section  5,  the 
one-dimensional  model  design  is  described  and  results  given.  Also 
in  this  section,  the  inverse  methodology  is  introduced,  applied, 
and  results  given.  Discussion  is  given  in  Section  6,  and  conclusions 
in  Section  7. 


2.  Model  description  and  verification 

2  1  SWAN  wave  prediction  model 

The  so-called  "third  generation”  (3G)  of  spectral  wave  models 
calculate  wave  spectra  without  a  priori  assumptions  regarding 
spectral  shape.  For  this  investigation,  we  use  the  SWAN  model 
("Simulating  WAves  Nearshore";  Booij  et  al.,  1999;  Holthuijsen 
et  al.,  2006).  SWAN  is  a  3G  model  designed  to  address  the 
excessive  computational  expense  of  applying  predecessor  3G 
models  (such  as  WAM,  WAMD1  Croup,  1988)  at  high  resolutions, 
particularly  in  coastal  regions.  The  governing  equation  of  SWAN 
and  most  other  3G  wave  models  is  the  action  balance  equation. 
In  Cartesian  coordinates,  the  action  balance  equation  is 

dN  dCxN  dCyN  dC„N  0C„N  S 

Of  dx  dy  Oa  d()  a 

where  a  is  the  angular  relative  frequency,  which  is  the  wave 
frequency  measured  from  the  frame  of  reference  moving  with 
current,  if  current  exists,  N  is  wave  action  density,  equal  to  energy 
density  divided  by  relative  frequency  ( N  =  E/a),  0  is  wave 
direction,  C  is  the  wave  action  propagation  speed  in  (x,  y,  a,  0) 
space,  e.g.  in  absence  of  currents,  Cx  is  the  x-component  of  the 
group  velocity  Cg ,  and  S  is  the  total  of  source/sink  terms  expressed 
as  wave  energy  density.  The  right-hand  side  of  the  governing 
equation  is  represented  by  three  terms,  S  =  Sm+Sr,,+Sds  (input  by 
wind,  nonlinear  interactions,  and  dissipation,  respectively).  The 
dissipation  term  can  be  broken  into  two  further  terms 
Sds  =  Sbr+Sf*,,;  the  Sbr  term  is  breaking  associated  with  steepness 
and  instability  (whitecapping,  surf  breaking,  etc.);  the  Sbo,  term 
includes  dissipation  due  to  bottom  roughness  Sbf  (friction, 
scattering),  percolation  S^,  or  non-rigid  bottoms  Smud  In  released 
versions  of  SWAN  (Holthuijsen  et  al.,  2006),  Sbot  is  only  associated 
with  rigid  seabeds,  S**,,  =  Sb/.  The  default  Sbf  formula  is  that 
of  JONSWAP  (Hasselmann  and  Coauthors,  1973),  in  which  the  user 
specifies  a  simple  tuning  coefficient  that  has  no  apparent  physical 
connection  with  measurable  seabed  characteristics.  An  alternate 
rigid-bed  formula  in  SWAN  is  that  of  Madsen  et  al.  (1988),  in 
which  the  user  specifies  a  single,  representative  bedform 
amplitude  at  each  point  in  the  computational  grid  on  which  Sbf 
is  estimated. 
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2.2.  Formulae  for  wove  damping  by  non-rigid  bottoms 

In  this  section,  three  Smud  formulae  are  described;  each  is  a 
different  representation  for  dissipation  by  a  viscous  fluid  mud 
layer.  All  three  formula  are  so-called  "two-layer"  models,  since 
they  represent  the  motion  of  both  the  mud  and  water  layers; 
however,  they  do  not  all  account  for  viscosity  in  both  layers.  These 
viscous  fluid  mud  models  vary  in  complexity  and  in  assumptions 
regarding  the  mud  thickness.  Analytical  models  which  account  for 
other  types  of  mud  behavior— viscoelastic  (e.g.  MacPherson,  1980; 
Jiang  and  Mehta,  1996;  Zhang  and  Ng,  2006)  or  plastic  (e.g.  Mei 
and  Liu,  1987) — are  not  included  in  this  study. 

2.2.1  Dalrymple  ond  Liu  ( 1978)  formulo 

Dalrymple  and  Liu  ( 1978)  (henceforth,  "DL")  propose  a  formula 
which  represents  wave  damping  due  to  viscosity  in  the  mud  layer 
and  the  overlying  water  layer.  It  is  a  relatively  accurate  approach, 
since  it  does  not  make  assumptions  about  the  thickness  of  the 
mud  layer.  It  requires  a  complex  iterative  solution  procedure  to 
compute  the  dissipation  rate.  DL  also  propose  a  modified  form, 
which  they  distinguish  from  the  other  as  a  "generally  applicable" 
technique  (in  their  Appendix  B),  using  the  assumption  that  the 
mud  layer  is  thin,  being  of  the  same  order  of  magnitude  as  the 
mud  Stokes’  boundary  layer;  this  thin-layer  model  also  requires 
iterations. 

2.2.2.  Ng  (2000)  formulo 

The  Ng  (2000)  (henceforth  "Ng")  formula  is  a  simplification  of 
the  DL  model.  There  are  two  major  differences  between  Ng  and  DL 
formulae.  Namely,  Ng  does  not  require  iterations  but  does  require 
that  the  mud  layer  is  much  thinner  than  the  overlying  water  layer. 
Like  DL  this  model  accounts  for  viscosity  in  both  layers. 

Ng  provides  an  expression  for  the  complex  wavenumber  k.  The 
imaginary  part  of  this  wavenumber,  k,  gives  the  wave  attenuation 
rate:  ?/(x,f)  =  Re[oe,(k*  0))  or  for  a  model  of  the  amplitude  decay, 

o  =  o,,e  k‘x  (1) 

where  t]  is  the  instantaneous  free  surface  elevation  and  o  is  the 
wave  amplitude.  We  implement  this  in  SWAN  by  considering  the 
case  of  a  single  wave  train  propagating  over  a  fiat  muddy  bottom, 
using  the  relation  (9/axXC**/V)/N  =  (S/ax)(Cg>xo2)/o2.  This  gives 
Sbot/E  =  -2Cgk,  (see  also  Komen  et  al.  (1994,  p.  170)).  Fig.  1  shows 


Fig.  1.  Variation  of  normalized  mud  sink  term  Sbo((fT,fl)/F(rT.tf)  (in  Hz)  according  to 
Ng  (2000)  as  implemented  herein.  Mud  thickness  <5mO  =  40cm  mud  density 
f>m  =  1310  kg/m3,  mud  kinematic  viscosity  vm  =  7.6  x  10  1  m2/s,  water  density 
()„  =  1000  kg/m3.  water  kinematic  viscosity  vw  —  1.0  x  10  6  m2/s. 


the  variation  of  Sbot/E  with  water  depth  and  wave  period  given 
values  of  density,  viscosity,  and  thickness  of  mud  employed  in 
Section  4  The  water  depth  range  (9-15  m)  also  corresponds  to  the 
range  for  which  mud  is  applied  in  Section  4. 

In  this  implementation  of  mud  effects  in  SWAN,  a  mud- 
adjusted  group  velocity  is  calculated  from  the  real  portion  of  the 
mud-adjusted  wavenumber.  Thus,  the  model  can  produce  “shoal- 
ing/de-shoaling”  effects  associated  with  spatial  variation  of  mud, 
separate  from  traditional  shoaling  associated  with  variation  of 
depths.  Mud-induced  "refraction"  has  also  been  added  but  this 
effect  is  not  included  in  the  present  study 

2.2.3.  Winterwerp  et  al.  (2007)  formula 

Another  formula  for  Smud  is  given  by  Winterwerp  et  al.  (2007) 
(henceforth  denoted  WDGL);  their  method  of  calculating  Smud  is 
based  on  Cade  (1958),  generalized  to  non-shallow  water  depths 
by  De  Wit  (1995).  Like  the  DL  method,  and  unlike  the  Ng  method, 
it  does  not  require  a  thin  mud  layer  and  does  require  iteration 
to  compute  the  dissipation  rate.  However,  the  iteration  in  the 
context  of  a  model  such  as  SWAN  is  not  computationally 
expensive.  Unlike  both  DL  method  and  Ng  method,  the  WDGL 
method  assumes  an  inviscid  water  layer.  However,  dissipation 
associated  with  water  viscosity  is  typically  quite  small 

In  WDGL  the  method  of  calculating  S^/E  from  k,  is  different 
from  the  implementation  of  Ng  used  in  this  study.  WDGL  follow 
Cade’s  method  (Equation  11-11  in  Gade  1958),  which  is  derived 
from  the  energy  transport  across  the  water/mud  interface, 
integrated  over  a  wave  period,  using  a  number  of  assumptions 
and  restrictions. 

In  the  version  of  the  Winterwerp  et  al.  (2007)  code  used 
herein,  the  effect  of  mud  on  phase  and  group  velocities  are  not 
considered.1 

2.4.  Verification  of  wave  damping  imp/emenfarions 

Fig.  2  compares  the  three  methods  for  estimating  the 
dissipation  rate  Dm  (Dm/W  in  the  case  of  DL  and  Ng).  which  is  the 
wave  attenuation  rate  due  to  viscous  mud  (plus  that  due  to 
viscous  water  in  the  case  of  DL  and  Ng),  equivalent  to  the 
imaginary  part  of  the  complex  wavenumber,  kr.  Variable  ks  here 
is  the  shallow  water  real  wavenumber,  rr/y/gh,  used  for  the 
normalization  consistent  with  DL  The  top  and  center  panels  use 
viscosity,  density,  water  depth,  wave  period  values  consistent  with 
Fig.  6  of  DL  The  DL  result  is  calculated  using  their  "thin  lower 
layer”  model.  The  lower  panel  uses  viscosity,  density,  water  depth, 
wave  period  values  consistent  with  Section  4  of  the  present  paper. 
Horizontal  axes  of  the  top  and  bottom  panels  are  the  normalized 
mud  layer  thickness,  d  =  dm$/Am,  where  <5m>0  is  the  fluidized  mud 
layer  thickness  and  Am  =-  y/2vm/(r  represents  the  Stokes  boundary 
layer  thickness. 

The  negative  impact  of  the  thin-layer  assumption  of  Ng  is 
noticeable  for  normalized  mud  depths  greater  than  approximately 
four.  For  smaller  values,  the  model  seems  quite  accurate.  In  both 
the  top  and  bottom  panels,  the  error  increases  with  larger  d ,  but  at 
a  different  rate,  suggesting  that  d  may  not  be  the  best  variable  to 
use  when  making  generalizations  regarding  the  limits  of  validity 
of  the  Ng  formula;  perhaps  Sm0  or  ^m.0/h  would  be  more 
appropriate. 

The  figure  shows  DL  thin  model  continuing  to  decrease  even 
for  large  values  of  d,  unlike  the  Ng  model.  The  top  panel  shows  the 
DL  thin  model  and  WDGL  diverging  for  large  d.  Unfortunately,  it  is 
not  known  which  is  more  accurate  in  this  range.  The  center  panel 


1  Subsequent  versions  do  have  this  feature  (Winterwerp.  personal  commu¬ 

nication  ). 
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Fig.  2.  Comparison  of  DL,  Ng.  and  WDC.L  (2007).  Top  panel:  comparison  with  DL  Fig.  6.  normalized  dissipation  as  a  function  of  normalized  mud  layer  thickness;  Center 
panel:  as  top  panel,  except  on  log  -log  scale,  and  variables  are  not  normalized.  Lower  panel:  as  top  panel,  but  with  different  settings,  chosen  to  be  more  representative  of  the 
Cassino  Beach  field  experiment. 


illustrates  the  effect  of  including  or  omitting  the  viscosity  of  the 
water  layer:  for  small  mud  thickness  values,  the  WDGL  model  for 
Dm  deviates  from  the  two  Dm/W  estimates,  but  this  occurs  only  in  a 
range  where  the  dissipation  is  very  small  by  any  method. 

In  Fig.  2,  top  and  bottom  panels,  there  exists  in  each  calculation 
method  a  value  of  J  for  which  the  dissipation  rate  is  a  maximum. 
There  is  negligible  difference  in  this  value  from  one  calculation 
method  to  another.  1.61,  1  57  and  1.57  for  DL  thin,  Ng  and  WDGL, 
respectively.  This  is  noteworthy  in  light  of  the  apparent  dis¬ 
crepancy  in  the  literature:  Gade  (1958)  suggests  a  value  of  1.2, 
while  Ng  gives  1.55.  The  text  of  DL  suggests  1.1  to  1.5,  but  this 
variation  is  apparently  due  to  simultaneously  using  two  defini¬ 
tions  for  ‘'peak  value  of  the  damping",  only  one  of  which  is 
consistent  with  our  definition. 

To  verify  the  implementation  of  Ng  and  WDGL  formulae  in  the 
SWAN  model,  the  expected  exponential  decay  from  (1)  is 
compared  with  the  actual  decay  calculated  by  SWAN.  For  the 
SWAN  simulations,  settings  are  identical  to  those  of  the  lower 
panel  of  Fig.  2  (and  consistent  with  Section  4),  and  all  energy  is 
contained  in  a  single  frequency/directional  bin.  This  is  shown  in 
Fig.  3.  The  results  show  that  the  k,  values  calculated  by  the  two 
methods  are  very  similar.  Further,  the  expected  decay  rate  with 
the  Ng  method  is  identical  to  the  Ng  SWAN  output.  However,  the 
WDGL  implementation  differs  from  its  corresponding  expected 
decay  rate.  This  is  due  to  the  different  method  of  calculatingSf,or/E 
from  k i  as  described  above.  In  separate  experiments,  it  was  found 
that  the  discrepancy  occurs  even  for  shallow  water  depths,  so  the 
shallow  water  assumption  of  Gade  (1958)  is  probably  not  the 
cause.  However,  the  methods  do  converge  for  large  wave  periods 
(e.g.  T  =  20  s).  We  do  not  assert  that  one  method  is  more  correct. 
However,  this  comparison  is  crucial  when  interpreting  the  results 
of  Section  4  and  5.  since  it  suggests  that  the  increased  dissipation 
with  WDGL  SWAN  versus  Ng  SWAN  in  the  field  application 


Fig.  3.  Comparison  of  decay  rates  for  Ng  (2000)  and  Winterwerp  et  al.  (2007). 

is  more  likely  due  to  differences  in  implementation,  rather  than 
differences  between  Ng  and  De  Wit  (1995). 

3.  Case  study  description 

3.1  Field  site  and  wave  climate 

The  field  site  of  Cassino  Beach  (Fig.  4)  was  chosen  due  to 
the  long-term  presence  of  a  large  offshore  mud  deposit  that 
periodically  transports  fluid  mud  throughout  the  nearshore  and 
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Fig.  4.  Map  of  the  study  area.  Fig.  4a,  Geographic  location  of  Cassino.  Fig.  4b. 
Cassino  region.  Outer  SWAN  grid  shown  in  red;  inner  SWAN  grid  shown  in  blue: 
transect  used  for  one-dimensional  model  shown  in  green. 


sometimes  onto  the  shoreline.  A  field  experiment  at  the  site 
occurred  in  May  and  June  2005.  The  bathymetry  in  the  vicinity 
of  Cassino  Beach  is  shown  in  Fig.  5,  with  locations  of  wave 
measuring  instruments  also  shown.  The  continental  shelf  slope 
here  is  very  mild,  approximately  25  m  vertical  over  35km 
horizontal  We  use  four  instruments  in  this  study  a  Datawell 
directional  waverider  buoy  at  approximately  25  m  water  depth,  a 
Nortek  acoustic  Doppler  profiler  (denoted  ”NDP")  at  8  m,  and  two 
pressure  gages:  "PA”  at  6  m  and  “P5”  at  2  m.  Description  of 
deployment  and  operation  of  the  pressure  gages  can  be  found  in 
Holland  et  al.,  this  issue. 

The  wave  climate  at  Cassino  Beach  is  dominated  by  windsea 
and  relatively  young  swells,  with  peak  period  of  8-12  s  being 
typical  Time  senes  of  significant  wave  height  from  the  previously 
mentioned  instruments  are  shown  in  Fig.  6.  Four  example 
directional  spectra,  derived  from  buoy  data,  are  also  shown  in 
the  figure,  A  large  majority  of  swell  fields  in  the  Southern  Ocean 
propagate  from  the  west  or  southwest.  Only  the  relatively 
uncommon  south-southwesterly  swells  from  the  Southern  Ocean 
propagate  in  the  direction  of  Cassino  Beach,  which  is  in  fact  most 
typically  sheltered  from  these  swells  due  to  the  concave  shape  of 
the  southeastern  South  American  coastline.  One  large  wave  event 


during  the  35-day  time  period  study  here  occurred  on  21-23  May 
2005  and  was  the  edge  of  a  fresh  swell  field  generated  by  a  large 
storm  south-southwest  of  Cassino  Beach,  east  of  Argentina  and 
Uruguay.  The  largest  wave  event— occurring  16-17  June  2005 — is 
from  a  small  but  intense  storm  just  offshore  of  Cassino  Beach. 

3.2.  Rheology 

Based  on  the  measurements  during  the  field  experiment 
(Holland  et  al.,  this  issue),  the  initial  best  estimates  of  rheology 
are  as  follows.  Note  that  even  though  the  experiment  data 
collection  capabilities  exceeded  prior  similar  attempts,  there  is 
considerable  uncertainty  in  these  estimates,  particularly  with 
respect  to  natural  deviations  from  measured  point  samples, 
therefore  an  additional  estimate  of  the  possible  range  of  values, 
is  given  in  parentheses: 

•  total  mud  layer  thickness,  =  40cm  (30- 100 cm), 

•  mud  kinematic  viscosity,  vm  =  7.6  x  10  3  m2/s  (1,4x10  3-15x 

10  3m2/s), 

•  mud  density.  pm  =  1310 kg/m3  (1080-1300  kg/m3), 

•  water  depth  where  mud  is  found:  h  =  9-15  m  (6  15  m). 

Here,  the  measured  total  mud  layer  thickness,  Sm.aj  is  given,  to 
distinguish  it  from  the  fluidized  mud  layer  thickness,  c5m0l  not 
available  from  the  observations  but  required  for  dissipation 
calculations. 

Prior  literature  has  drawn  attention  to  very  fast  wave 
attenuation  possible  due  to  viscous  mud.  However,  this  rate 
depends  heavily  on  mud  characteristics  and  water  depth.  To 
demonstrate  this,  Table  1  compares  dissipation  rates  and  relevant 
variables  for  five  scenarios,  including  two  previously  published 
works,  Gade  (1958)  and  Dalrymple  and  Liu  (1978).  For  the  latter, 
the  example  used  in  DL’s  discussion  of  their  Fig.  8  is  used  here. 
Kaihatu  et  al.  (2007)  refers  to  a  manuscript  in  which  the  Ng 
formulation  is  applied  in  a  phase-resolved  model.  “Cassino  (2005 
example)”  uses  the  environmental  characteristics  taken  from 
observations  from  the  field  experiment  studied  in  the  present 
paper,  assuming  that  the  entire  mud  layer  is  fluidized,  i.e. 
^m.o  =  A  hypothetical  Cassino  Beach  scenario  is  presented, 
identical  to  the  2005  case  except  that  the  mud  is  located  at  or  near 
the  surf  zone  The  exponential  decay  rate  is  illustrated,  here,  by 
showing  the  distance  of  propagation  at  which  a  wave  would 
be  attenuated  to  50%  or  10%  of  its  original  amplitude,  indicated 
as  x<j  osao  and  xfl=oia0  in  the  table;  this  calculation  assumes  a 
flat  bottom,  with  no  other  source/sink  terms  active  (e.g.  a  low 
steepness  wave  with  no  wind). 

It  is  immediately  apparent  that  the  hypothetical  field  cases 
used  by  DL  and  Kaihatu  et  al.  (2007)  shows  very  fast  attenuation 
relative  to  either  Cassino  example.  All  the  three  non-Cassino 
examples  have  d  between  1.1  and  1.6,  a  range  in  which  the 
dissipation  by  mud  is  expected  to  reach  a  maximum  value  relative 
to  thinner  or  thicker  mud  layers.  In  the  Cassino  examples, 
d  =  2  57,  well  above  this  peak;  for  this  value  of  relative  thickness, 
slightly  decreasing  the  mud  thickness  would  be  expected  to 
increase  dissipation  rate;  a  large  decrease  in  mud  thickness,  say 
to  <5m0  =  5  cm,  or  any  increase  in  mud  thickness  is  expected  to 
decrease  the  dissipation  rate.  Also  noteworthy  is  the  high  viscosity 
used  by  DL,  500  times  greater  than  that  used  for  our  Cassino 
application.  As  expected,  in  the  Cassino  example,  when  the  mud 
is  at  or  near  the  surf  zone,  the  dissipation  is  much  faster:  xa  05ao  = 
330m  compared  to  5-6  km  with  the  mud  observed  further 
offshore,  as  it  was  during  the  2005  field  experiment. 

For  the  2005  Cassino  example,  sensitivities  of  dissipation  level 
to  wave  period  and  water  depth  are  illustrated  previously  in  Fig.  1. 
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Fig.  5.  8dthymetry  and  mud  layer  applied  in  two-dimensional  simulations  Left  panel:  bathymetry  on  outer  SWAN  grid,  with  nest  location  and  buoy  location  indicated 
Center  panel:  location  of  applied  mud  layer  in  outer  SWAN  grid.  Upper  right  panel:  bathymetry  on  inner  SWAN  grid  (nest),  with  horizontal  position  of  wave-measuring 
instruments  shown  Lower  right  panel*  depth  profile,  from  origin  to  buoy,  with  P5,  PA,  NDPand  buoy  locations  indicated 
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Fig.  6.  Time  series  of  significant  wave  height  observed  at  four  locations.  "WR"  denotes  the  waverider  buoy.  Numbers  in  parentheses  are  approximate  local  depths  in  meters. 
Directional  wave  spectra,  inferred  from  buoy  data,  are  also  shown,  corresponding  to  four  time  periods:  0007  UTC  22  May,  0337  UTC  24  May.  1139  UTC  29  May.  and  1209  UTC 
16  June. 


Table  2  illustrates  the  sensitivity  to  variation  in  mud  viscosity, 
density  and  thickness,  using  the  possible  range  of  values  given 
above.  Also  included  is  a  calculation  for<5m0  =  10cm,  allowing  for 
the  possibility  that  <5m.o><Vo.r-  It  is  apparent  that  for  this  wave 
period  and  water  depth  (Tp  =  10s  and  h  —  12  m),  uncertainty  in 
the  viscosity  is  of  concern,  with  decay  length  scale  varying  from 
Xq  o  5ao  ~  3.6  to  12.5  km  for  the  highest  and  lowest  viscosity 
values,  respectively.  Further,  the  computations  indicate  that  if 


only  a  fraction  of  the  mud  layer  is  fluidized,  this  will  also  have  a 
large  impact.  Sensitivity  to  mud  density  is  relatively  minor. 

4.  Two-dimensional  modeling 

In  this  section,  we  apply  SWAN  model  with  and  without 
dissipation  by  viscous  mud.  This  is  performed  without  any  tuning 
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Table  1 

Comparison  of  dissipation  rates  and  relevant  variables  for  five  scenarios 


Cade  (1958)  Fig.  4 

Datrymple  and  Liu  (1978)  Fig.  8 

Xaihatu  et  at.  (2007) 

Cassino  (2005  example) 

Cassino  (hypo-thetical) 

Tp(  s) 

1.40 

5.0 

100 

10.0 

10.0 

h  (m) 

0,038 

4.0 

1.0 

12.0 

2.0 

pm  (kg/m3) 

1504.00 

1800.00 

1111.11 

1310.00 

1310.00 

y 

0  57 

0.57 

0  90 

0.76 

0.76 

vw  (m2/s) 

2.42e-06 

2.60e-06 

1.00e-06 

1.00e-06 

I.OOe-06 

vm  (m2/s) 

2.60e-03 

1.49 

1.00e-02 

7.60e-03 

7,60e-03 

C 

32.78 

756.26 

100  00 

87.18 

8718 

&mjo  ( m) 

0.038 

2.0 

0.20 

0400 

0.400 

•Mm) 

0.03 

1.54 

018 

016 

0.16 

2 

1.12 

1.30 

1.12 

2.57 

2  57 

Dm/w  =  k,  (rad/m) 

0.90 

9.81e-03 

7.51  e-03 

1.29e-04 

2.09e-03 

Vo soo  (m) 

0.77 

70.65 

92.29 

5364 

332 

Xo-o too  (m) 

2.56 

234.7 

306.6 

17821 

1103 

All  dissipation  rates  here  are  computed  using  Ng  (2000). 


Table  2 

Sensitivity  to  mud  viscosity,  density  and  thickness,  given  Tp  =  10s  and  h  =  12  m, 
according  to  Ng  (2000) 


(m2/s) 

Pm  (kg/m3) 

<*m.0  (m) 

Dmiw~  M rad/m) 

Xo^osoo  (km) 

1  4e-3 

1310 

0.40 

5.50e-5 

12.5 

7.6e-3 

1.29e-4 

5.4 

15e-3 

1.94e-4 

3.6 

7.6e-3 

1080 

040 

1.57e-4 

4.4 

1310 

1.29e~4 

5.4 

7  6e-3 

1310 

0.10 

4.34e-S 

16.0 

0  30 

1.36e-4 

SI 

0.40 

1.29e-4 

5.4 

1.00 

1.29e-4 

5.4 

to  improve  agreement  at  the  three  onshore  measurement 
locations.  A  total  of  431  simulations  are  used  to  evaluate  trends 
and  sensitivity  to  forcing  conditions,  and  to  generate  statistics. 
The  simulations  without  dissipation  by  viscous  mud  help  address 
the  objective  of  evaluating  the  necessity  of  (or  lack  thereof) 
including  physics  of  wave  dissipation  by  viscous  mud  in  this 
application.  The  simulations  with  dissipation  by  viscous  mud  are 
performed  using  the  Ng  SWAN  and  WDGL  SWAN  models;  these 
simulations  address  another  objective,  which  is  to,  given  our  best 
guess  of  rheology  and  mud  distribution  derived  from  field 
measurements,  evaluate  the  skill  of  these  models  to  predict 
observed  wave  heights. 

Directional  spectra  for  boundary  forcing  were  available  for  746 
time  periods  between  May  15  and  June  25.  The  process  for 
reducing  this  population  to  431  time  periods  is  summarized  now 
14  of  746  boundary  spectra  were  eliminated  in  quality  checking 
(Section  4.1);  163  of  the  remaining  732  were  discarded  due  to 
absence  of  reliable,  coincident  wind  data  (Section  4.2);  46  of  the 
remaining  569  were  discarded  due  to  absence  of  coincident  NDP 
data;  528  time  periods  were  simulated;  97  of  these  528  were 
discarded  due  to  mismatch  at  the  offshore  boundary  (Section  4.6), 
leaving  431  time  periods  for  calculation  of  statistics. 

4.1  Boundary  forcing 

The  wave  model's  outer  computational  grid  is  designed  such 
that  the  center  of  the  offshore  boundary  corresponds  to  the 
location  of  the  Datawell  Waverider  buoy,  from  which  directional 
measurements  are  available  for  the  time  period  of  interest. 

Directional  data  were  provided  in  the  form  of  five  variables 
given  at  each  of  64  frequencies  non-uniformly  spaced  from  0.025 
to  0.58  Hz.  The  five  variables  are  spectral  energy  density  and  four 
directional  moments  describing  the  unmeasured  normalized 


directional  distribution  for  that  frequency.  D(0).  This  information 
was  available  for  746  time  periods  between  May  15  and  June  25. 
The  procedure  for  estimating  D(0)  is  as  follows:  The  four  Fourier 
coefficients  corresponding  to  these  four  moments  were  calculated 
directly.  Then,  an  estimate  of  D(0)  was  calculated  using  the 
maximum  entropy  method  (Lygre  and  Krogstad,  1986).  In  the  case 
of  14  spectra,  the  integrated  spectra  did  not  match  the  waveheight 
in  the  header  of  the  data  file;  these  were  discarded.  Thus,  732 
E((J,0)  estimates  were  retained  and  converted  to  a  format  readable 
by  SWAN 

4.2.  Wind  /orcing 

Wind  vectors  were  estimated  using  high-frequency  energy 
measured  by  the  waverider  buoy.  For  each  of  732  spectra,  the 
wind  speed  is  chosen  that,  in  fetch-limited  conditions,  results  in  a 
SWAN-predicted  high-frequency  spectrum  consistent  with  the 
measured  spectrum.  This  algorithm  was  able  to  create  a  relatively 
confident  prediction  for  569  of  732  cases.  The  algorithm  is 
described  in  Appendix  A. 

4.3.  Bathymetry  and  wafer  level  specification 

Regional  charts,  ship  surveys,  digitized  shorelines,  and  near¬ 
shore  bathymetry  from  a  jet-ski  system  were  used  to  compute  the 
bathymetric  surface  used  in  the  model.  The  data  were  corrected 
to  remove  tidal  offsets  and  merged  using  the  method  of  Plant  et  al. 
(2002).  The  nearshore  surveys  were  made  during  late  April  and 
May  2005.  See  Holland  et  al  (this  issue)  for  further  description. 

The  tide  range  at  Cassino  Beach  is  approximately  ±60 cm  from 
the  mean  water  level,  with  a  significant  meteorological  compo¬ 
nent.  Bathymetry  is  relative  to  mean  lower  low  water;  thus,  an 
offset  is  determined  for  each  hindcast  using  pressure  gage  data. 
For  the  period  26  May-24  June,  "P2"  in  approximately  1.2  m  water 
depth,  corrected  for  barometric  pressure,  is  used.  Prior  to  26  May, 
the  "PA”  gage  is  used,  also  used  in  the  wave  comparisons  below. 

4.4  Mud  specification 

For  the  two-dimensional  application  of  the  Ng  SWAN,  we  used 
the  mud-related  variables  given  of  the  initial  best  estimate 
[<5m  0  =  40cm,  vm  =  7.6e-03  m2/s,  pm  =  1310  kg/m3,  as  described 
in  Section  3.2).  For  the  horizontal  location,  the  mud  is  applied 
where  the  water  depth  is  between  9  and  15  m  (see  center  panel  in 
Fig.  5).  In  locations  where  mud  is  not  applied,  the  source  term  S 
is  assumed  zero.  In  other  words,  bottom  friction  associated  with 
sand  or  otherwise  rigid  seafloor  is  not  included.  This  is  motivated 
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by  lack  of  measurements  describing  possible  rigid  bedforms. 
(Separate  simulations  with  dissipation  by  rigid  bed  forms  via 
Madsen  et  al.  (1988)  were  produced;  these  results  are  not 
presented  here,  but  are  discussed  in  Section  6.) 

4.5.  Grid  specification,  wave  generation  physics 

An  inner  nest  was  included  in  the  simulations  to  provide 
higher  resolution  in  the  nearshore;  the  nest  position  is  shown  in 
Figs.  4b  and  5.  Additional  details  on  the  model  design  are  as 
follows: 

•  a  rotated  rectilinear  grid,  such  that  the  positive  x-direction  is 
48.33  clockwise  from  the  west-to-east  direction, 

•  for  the  x-axis,  shore-normal  and  positive  in  the  southeast 
direction,  the  grid  spacing  Ax  =  252  m,  with  155  grid  nodes 
(and  for  the  inner  nest,  50  m  and  146  nodes), 

•  for  the  y-axis,  shore-parallel  and  positive  in  the  northeast 
direction,  the  grid  spacing  Ay  —  500  m,  with  402  grid  nodes 
(and  for  the  inner  nest,  50  m  and  365  nodes), 

•  in  frequency  space,  34  bins  in  logarithmic  distribution  from 
0.0418  to  1.0  Hz, 

•  10  directional  resolution  (36  directional  nodes), 

•  stationary  computations,  with  default  settings  for  numerics, 

•  default  settings  for  nonlinear  interactions  Sn(4  and  depth- 
limited  breaking  Sdsbn 

•  for  whitecapping  term,  Sds<vvo  the  van  der  Westhuysen  et  al. 
(2007)  formula  is  used  to  reduce  problems  associated  with 
non-physical  dependence  on  mean  steepness  (Rogers  et  al., 
2003  and  references  therein),  and 

•  for  wind-to-wave  energy  transfer  term,  Sin,  the  Yan  (1987) 
formula  is  used  (see  van  der  Westhuysen  et  al.  (2007)). 


4.6.  Computations  and  Results 

Five  hundred  and  twenty  eight  two-dimensional  simulations 
with  each  of  the  three  SWAN  variants  were  performed  (Shot  =  0, 
Shot  =  Smud  according  to  Ng.  and  S^t  =  Smud  according  to  WDGL). 
Simulations  were  discarded  if  any  of  the  three  models  produced  a 
mismatch  of  wave  height  at  the  offshore  boundary  (greater  than 
4%  normalized  error),  indicative  of  problems  with  wave  growth 
internally.* 2  Ninety  seven  time  periods  were  omitted,  leaving  431 
cases.  Results  for  two-dimensional  simulations  are  shown  in  Fig.  7 
and  Table  3.  Statistics  shown  in  the  table  are:  number  of 
comparisons,  bottom-induced  dissipation  applied,  bias,  rms  error, 
correlation  coefficient,  standard  deviation  of  error  (i.e.  rms  error 
with  bias  effect  removed),  scatter  index,  and  the  mean  of  observed 
values.  The  scatter  index  is  defined  as  (e.g.  Cardone  et  al.  1995) 


and  correlation  coefficient  is 

cc  (0-6)(M-M) 
v/(O-0)2v/Tm-M)2 

where  overscore  indicates  a  mean.  0  are  observations  and  M  are 
model  values. 

It  is  apparent  in  the  comparisons  for  P5  that  a  distinct  sub- 
population  exists  for  which  observed  waveheights  are  very  small 

2  This  is  a  surprising  result,  since  the  wind  speed  inference  algorithm  should 
provide  the  optimal  wind  speed  for  producing  in  the  model  the  same  wind  sea  as 
observed  at  the  buoy.  This  is  due  to  non-physical  sea-swell  interaction  and  is 
discussed  further  in  the  Appendix  A. 


(less  than  38  cm)  and  much  smaller  than  waveheights  from  either 
of  the  two  models.  These  correspond  to  two  intervals  during  June 
6-15  It  was  determined  that  the  small  waveheights  are  very  likely 
due  to  instrument  failure  (the  gage  may  have  been  covered  by 
sediment)  and  therefore  these  points  are  not  included  in  the 
calculated  statistics. 

For  the  simulations  without  mud-induced  dissipation  (left 
column  in  Fig.  7),  there  is  a  clear  trend  of  overprediction  of  wave 
energy,  with  bias  of  33,  46,  and  16  cm  at  NDP,  PA  and  P5, 
respectively.  This  addresses  one  of  the  previously  stated  objec¬ 
tives,  insofar  as  it  suggests  that  some  type  of  bottom-induced 
dissipation  is  needed  for  accurate  prediction  of  the  nearshore 
energy  level.  Implied  in  this  conclusion  is  the  assumption  that 
there  is  no  other  source  of  error — in  the  observations  or 
simulations — that  would  lead  to  large  positive  bias;  possible 
errors  are  discussed  in  Section  6.  A  surprisingly  small  bias  exists 
for  the  high-energy  simulations  without  mud-induced  dissipa¬ 
tion.  In  fact,  there  is  little  or  no  bias  for  the  time  periods  with 
largest  waveheight  in  the  NDP  and  P5  records.  One  possible 
explanation  is  that  the  viscosity  could  be  reduced  under  stronger 
level  of  forcing,  as  might  be  expected  with  a  thixotropic  (shear¬ 
thinning)  fluid  Or,  the  water/mud  interface  may  be  indistinct 
during  high  energy  events,  e.g.  due  to  sediment  suspension, 
preventing  effective  energy  transfer. 

In  the  simulations  with  mud-induced  dissipation  (right 
column  in  Fig.  7),  the  positive  bias  is  eliminated.  At  the  two 
shallower  locations,  PA  and  P5,  the  bias  with  the  Ng  model  is  very 
small  relative  to  the  S*,,  =  0  simulations  (-6  and  -8  cm, 
compared  to  +46  and  +16  cm).  At  the  NDP  location,  the  bias  is 
similar,  but  of  opposite  sign  (-29cm  compared  to  +33  cm). 
The  negative  bias  at  this  location  suggests  that  the  modeled 
dissipation  is  too  strong  in  the  region  between  the  NDP  and  the 
waverider  buoy,  with  the  implied  assumption  as  noted  above.  One 
explanation  is  that  the  true  mud  distribution  in  this  region  might 
be  more  thin,  patchy,  or  less  fluidized  than  the  uniform  lens  of 
40cm  fluid  mud  applied  here;  this  is  explored  further  via  inverse 
modeling  in  Section  5.  The  rms  error  and  scatter  index  are  mostly 
improved  by  including  Smud,  while  the  correlation  is  not  always 
improved;  in  fact,  it  is  significantly  worse  at  PA. 

Comparing  the  Ng  and  WDGL  statistics,  it  is  clear  that  the 
dissipation  is  generally  stronger  with  the  latter.  Considering  the 
comparisons  in  Section  2,  this  is  most  likely  due  to  difference  in 
implementation  rather  than  in  the  formulae  themselves. 


5.  One-dimensional  modeling:  forward  and  inverse  solutions 

In  the  above-mentioned  two-dimensional  modeling,  the  model 
with  Shot  =  0  tends  to  overpredict  wave  energy  and  the  models 
that  include  Smud  tend  to  underpredict  energy.  Assuming  that 
modeling  errors  unrelated  to  Sbot  are  relatively  small,  this  result 
suggests  that  the  wave  dissipation  Smud  occurred  at  Cassino  but  is 
overpredicted  by  the  model  if  the  calculations  are  made  using 
mud  characteristics  based  on  field  measurements  (Section  3.2). 
The  range  of  uncertainty  in  the  measured  mud  viscosity  and  the 
thickness  of  the  fluidized  mud  layer  (Table  2)  are  not  small,  and 
the  mud  distribution  shown  in  Fig.  5  is  unlikely  to  be  an  accurate 
portrayal  of  the  heterogeneity  of  the  actual  mud  deposits.  Further, 
the  mud  distribution,  thickness,  and  depth  of  fluidization  may  not 
have  been  constant  during  the  simulation  period,  core  samples — 
which  in  fact  only  identified  the  first  two  of  these  three 
variables— were  taken  during  only  one  week  of  May  2005.  The 
dissipation  calculations  given  in  Table  2  suggest  that  the  model 
results  could  be  sensitive  to  these  uncertainties  in  the  forcing. 

One  method  of  quantifying  the  sensitivity  of  the  results  to  the 
uncertainty  in  forcing  is  to  run  the  forward  model  with  a  number 
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Fig.  7.  Scatter  plot  comparison  of  SWAN  vs.  observations  at  three  locations,  both  with  and  without  dissipation  by  viscous  mud  included  in  SWAN  (Ng  formula). 


of  forcing  conditions  and  compare  the  results.  However,  due  to  the 
non-monotonic  behavior  evident  in  Fig.  2,  a  large  number  of 
hypothetical  conditions  would  be  required  to  quantify  the 
sensitivity.  A  more  deterministic  approach  is  preferred.  To  this 
end,  inverse  modeling  was  used  to  determine  if  any  mud 
distribution  would  yield  the  observed  waveheights,  and  if  so, 
what  that  distribution  is.  (Sensitivity  to  uncertainty  in  viscosity 
was  also  explored  in  this  manner;  these  results  are  discussed 
qualitatively  in  Section  6.) 

The  forward  modeling  was  a  blindfold  process,  so  the 
NDP,  PA,  and  P5  wave  observations  were  disregarded 
until  the  comparison  stage.  Similarly,  the  inversion  for  the  mud 
distribution  does  not  utilize  much  of  the  relevant  knowledge 
obtained  from  the  mud  observations.  Most  importantly,  the 
inversion  for  mud  distribution  permits  solutions  suggesting  mud 


at  locations  shallower  than  the  9  m  depth  contour,  which  is 
not  consistent  with  existing  observations.  Correspondence, 
or  lack  thereof,  between  the  inversion  results  and  the  mud 
observations  provides  insight  unavailable  from  the  forward 
modeling  alone. 

Inverse  modeling  requires  multiple  simulations  for  each  of  the 
time  periods,  so  creating  inverse  solutions  for  each  of  431  time 
periods  might  normally  require  more  computation  time  than  is 
practical.  This  challenge  was  addressed  by  reducing  the  two- 
dimensional  simulations  to  one  geographic  dimension,  specifi¬ 
cally  a  cross-shore  transect  running  from  the  waverider  buoy 
to  the  shoreline.  Of  course,  by  ignoring  long-shore  variation  in  the 
system,  some  error  is  introduced,  so  in  Section  5.2  below,  the 
forward  two-dimensional  model  is  compared  with  the  forward 
one-dimensional  model. 
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Table  3 

Significant  wave  height  statistics  for  two-dimensional  simulations 


n 

s** 

8ias 

(m) 

rmse 

(m) 

CC 

side 

(m) 

SI 

<obs> 

(m) 

WR 

431 

0 

0.01 

0.04 

1  00 

0.03 

0.02 

1  72 

Ng  (2000) 

0.01 

0.03 

1.00 

0.03 

002 

- 

WDGL 

001 

0.03 

1.00 

003 

002 

NDP 

431 

0 

033 

0.37 

097 

0.17 

0.36 

1.03 

Ng  (2000) 

-0.29 

0.38 

0.96 

0.2S 

0.37 

- 

WDGL 

-0.40 

0.S1 

0.94 

0.31 

0.49 

PA 

55 

0 

0  46 

0.48 

0.94 

0.13 

074 

0.65 

Ng  (2000) 

-0.06 

0.17 

79 

0.17 

0.27 

- 

WDGL 

-0.12 

0.20 

0.81 

016 

0.31 

P5 

323 

0 

0  16 

022 

0  84 

0.15 

0.29 

0.76 

Ng (2000) 

-0,08 

0.11 

0.95 

0.09 

0,15 

- 

WDGL 

-0.1S 

0.18 

0.93 

0.09 

0.24 

- 

5.1  One-dimensional  grid  specification 

Model  design  choices  differing  from  the  two-dimensional 
simulations  are  as  follows: 

•  For  the  x-axis,  shore-normal  and  positive  in  the  southeast 
direction,  the  grid  spacing  /lx  =  50  m.  with  771  grid  nodes 
(no  inner  nest). 

•  There  is  no  shore-parallel  y-axis. 

•  A  curvature-based  stopping  criterion  is  used,  thus  typically 
computing  more  iterations  than  would  be  the  case  if  the 
default  criterion  was  used.  This  potentially  improves  accuracy, 
see  Zijlema  and  van  der  Westhuysen  (2005). 


5.2.  Forward  model  results 

Table  4  compares  the  one-dimensional  forward  model  results 
for  several  models: 

•  no  dissipation  by  mud, 

•  dissipation  by  mud  with  40cm  mud  layer  thickness  applied 
between  9  and  15  m  water  depth  contours  using  Ng  SWAN, 

•  dissipation  by  mud  with  40cm  mud  layer  thickness  applied 
between  9  and  15  m  water  depth  contours  using  WDGL  SWAN, 

•  dissipation  by  mud  with  1  m  mud  layer  thickness  applied 
between  9  and  15  m  water  depth  contours  using  WDGL  SWAN, 
and 

•  inversion  results,  with  mud  thickness  adjusted  to  retrieve 
observed  wave  heights. 

Statistics  for  the  two-dimensional  and  one-dimensional  simu¬ 
lations  can  be  compared  in  Tables  3  and  4.  It  is  apparent  that  there 
are  some  differences  between  the  one-  and  two-dimensional 
results,  but  none  that  would  actually  change  the  conclusions.  The 
most  noteworthy  difference  is  in  the  bias  of  Ng  SWAN  at  PA,  which 
changes  from  small  negative  bias  in  the  two-dimensional  case 
(-6cm)  to  slightly  positive  in  the  one-dimensional  case  (+1  cm). 

Knowing  that  application  with  r)m0  =  40cm  results  in  too 
much  dissipation  (i.e.  positive  wave  height  bias),  it  is  useful  to 
consider  whether  application  with  <5m0  =  1  m  will  improve  this 
bias.  For  a  typical  boundary  layer  thickness  value,  say  Am  =  16  cm 
(Table  1),  the  mud  thickness  of  ^m<0  =  40cm>  1.6. 1m.  Referring  to 
Fig.  2,  this  suggests  that  increasing  the  mud  thickness  will  result 
in  a  decrease  in  dissipation  and,  therefore,  improve  the  negative 


Table  4 

Significant  wave  height  statistics  for  one-dimensional  simulations 


n 

8ias 

(m) 

rmse 

(m) 

CC 

side 

(m) 

SI 

<obs> 

(m) 

WR 

431 

0 

002 

0.05 

1.00 

0.05 

0.03 

1.72 

Ng  (40cm) 

0.02 

O.OS 

1.00 

0.04 

0.03 

- 

WDGL  (40cm) 

0.02 

0,05 

1.00 

0.04 

0.03 

- 

WDGL  (1  m) 

0.02 

O.OS 

1.00 

0.04 

0.03 

- 

WDGL  Inverse 

0.02 

0.05 

1.00 

0.04 

0.03 

NDP 

431 

0 

036 

0.41 

096 

0.19 

0.39 

1.03 

Ng (40cm) 

-0.22 

0.31 

0.96 

0.22 

0.30 

- 

WDGL  (40  cm) 

-0.32 

0.42 

0.94 

0.28 

0.41 

- 

WDGL  (1  m) 

-0.31 

0.41 

0.94 

0.27 

0.39 

- 

WDGL  Inverse 

-0.02 

0.03 

1.00 

0.02 

0.03 

PA 

SS 

0 

0.52 

0.55 

092 

0.18 

0.8S 

0.6S 

Ng  (40cm) 

0.01 

0.16 

0.77 

0.17 

0.25 

- 

WDGL  (40cm) 

-0.04 

0.17 

0.79 

0.16 

0.26 

- 

WDGL  ( 1  m) 

-0.03 

0.16 

0  80 

0.16 

0.24 

- 

WDGL  Inverse 

0.01 

0.03 

1.00 

0.02 

0.04 

P5 

323 

0 

0.15 

0.22 

0.79 

0.16 

0.29 

0.76 

Ng  (40cm) 

-0.06 

0.10 

0.9S 

0.08 

0.13 

- 

WDGL  (40cm) 

-0.12 

0.14 

0.95 

0.08 

0.18 

- 

WDGL  ( 1  m) 

-0.11 

0.13 

0.95 

0.08 

0.17 

- 

WDGL  Inverse 

-0.02 

0.04 

0.99 

0.04 

0.06 

- 

bias  seen  for  the  40  cm  WDGL  set  in  Table  4.  However,  the  bias  for 
the  1  m  WDGL  set  is  only  slightly  better  (for  NDP,  a  change  from 
-29  to  -28cm),  indicating  that  the  dissipation  in  this  thickness 
range  is  only  weakly  sensitive  to  mud  thickness  (consistent  with 
Fig.  2),  and  zero  bias  cannot  be  achieved  without  using  a  fluidized 
mud  layer  thickness  <5,n0  which  is  less  than  the  probable  range  of 
total  mud  thickness  (<5m>0 j  =  30cm-l  m,  Section  3.2)  or  a  smaller 
mud  viscosity. 

5.3.  Inverse  modeling  method  description 

Given  the  uncertainties  in  the  forward  modeling  outlined 
above,  an  inversion  methodology  was  developed.  The  inversion  is 
simple:  the  shore-normal  transect  which  constitutes  the  compu¬ 
tational  grid  is  split  into  zones,  and  working  in  a  shoreward 
direction,  the  mud  layer  thickness  is  determined  which  results  in 
the  wave  height  that  was  observed  at  the  shoreward  end  of  each 
zone.  The  zones  are  as  follows: 

( 1 )  Zone  A,  from  the  offshore  boundary  to  15  m  water  depth.  It  is 
assumed  that  no  mud  exists  here,  <5m0  =  0. 

(2)  Zone  B,  from  15  m  water  depth  to  7.3  m  depth,  terminating  at 
NDP 

(3)  Zone  C,  from  7.3  to  5.6  m  depth,  terminating  at  PA. 

(4)  Zone  D,  from  5.6  to  1.8  m  depth,  terminating  at  P5. 

(5)  Zone  E,  from  1.8  m  depth  to  the  shoreline:  Sm0  =  0  assumed, 
but  dissipation  level  in  this  zone  is  actually  irrelevant  to  the 
tabulated  model-data  comparisons,  since  waves  are  not 
measured  after  they  pass  through  this  zone. 

In  cases  where  PA  data  are  not  available,  Zones  C  and  D  are 
combined  and  inversion  for  the  two  zones  is  based  on  waveheight 
at  P5.  In  cases  where  P5  data  are  not  available,  Zone  D  is  treated 
similar  to  Zone  E,  i.e.  <5m,o  =  0  assumed. 

One  inherent  assumption  in  the  inversion  is  that  the 
waveheight  as  an  observation  point  is  affected  only  by  the  mud 
offshore  of  that  point,  e.g.  the  waveheight  at  NDP  is  not  affected 
by  the  mud.  or  lack  thereof,  in  Zone  C.  This  assumption  allows  the 
procedure  to  advance  in  a  stepwise  fashion,  from  Zone  B  to  C  to  D, 
rather  than  solving  all  three  simultaneously.  Since  the  slope  is 
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very  mild,  and  there  are  many  cases  with  windsea  propagating 
toward  deeper  water  (i.e.  offshore  winds),  this  assumption  is  not 
obviously  valid.  The  underlying  physical  assumption  is  that  these 
fetch-limited  waves  are  too  short  to  induce  significant  pressures 
near  the  bottom  and,  therefore,  induce  motion  in  the  mud  layer. 
The  validity  of  this  assumption  is  discussed  further  in  Section  5.4. 

Based  on  the  similarity  of  the  two  models  in  the  previous 
results,  there  was  no  compelling  reason  to  perform  the  inversion 
for  both  Ng  SWAN  and  WDCL  SWAN.  The  latter  was  chosen  for  the 
inversion,  since  that  model  is  valid  for  a  larger  range  of  mud 
thickness  (Fig.  2). 

An  inverse  solution  is  possible  when  the  observed  wave  height 
is  bounded  by  the  result  with  zero  mud  and  the  result  with  most 
dissipative  possible  mud  thickness:  <Hobs 

The  thickness  <5m.0.mjx  typically  corresponds  to  15<<?<16. 
i.e.  T5/Im  in  our  applications  of  WDGL  and  Ng. 

The  relationship  of  dissipation  (and  therefore,  waveheight  in  the 
inversion)  with  mud  thickness  always  has  a  single  peak  (trough). 
Thus,  when  one  solution  exists,  there  is  always  one  other  solution. 
This  is  illustrated  in  Fig.  8,  where  the  horizontal  line  represents 
the  observed  waveheight  that  the  inversion  should  produce.  Since 
the  lines  intersect,  there  exists  a  pair  of  solutions:  one  solution 
exists  at  7  cm  and  another  at  11  rn,  It  is  apparent  that 

in  the  thick  region  (which  is  the  right-hand  side  of  the  figure, 
corresponding  to  <5m 0<  1.5zfm),  the  dissipation,  and  therefore  the 
resulting  waveheight.  is  only  weakly  sensitive  to  the  mud 
thickness.  When  two  solutions  exist,  the  procedure  must  choose 
one.  In  the  present  implementation,  the  inversion  always  chooses 
the  thinner  solution,  since  the  thicker  solution  was  often 
unrealistic,  as  in  this  example  (<5m0  =  11  m).  The  waveheight  from 
the  forward  model  is  not  a  linear  function  of  mud  thickness;  thus, 
the  problem  must  be  solved  via  an  iterative  sequence  of  forward 
simulations.  The  solid  line  in  Fig.  8  shows  the  waveheight  from 
many  (119)  applications  of  the  forward  model,  each  using  a 
different  mud  thickness.  In  practice,  this  many  forward  model 
runs  is  not  necessary;  rather,  a  modified  Newton  Raphson 
procedure  is  used  to  find  the  solution  within  a  relatively  small 
number  of  forward  simulations.  When  the  problem  is  not 
bounded  (i.e.  the  inversion  cannot  match  the  observed  wave- 
heights  exactly,  which  would  be  the  case  if  the  two  lines  in  Fig.  8 
did  not  intersect)  the  solution  is  used  which  yields  the 
smallest  waveheight  mismatch;  specifically  when  > Hobs, 

dm, 0  =  ^m.O.max  ^tld  When  Hobs  >H^m0= q,  ^rn.O  0- 


mud  layer  thickness,  <5m0  (m) 


Fig.  8  Example  of  inferring  mud  thickness  from  wave  height  observation.  Here, 
solutions  of  <Vo  =  0  07  and  11  m  are  obtained.  In  this  case,  the  maximum 
dissipation  occurs  with  18.4cm. 


5.4.  Inverse  model  results 

The  inversion  solution  yielded  an  exact  match  in  waveheight  in 
all  zones  for  over  half  (248)  of  the  431  cases.  In  the  other  cases,  the 
error-minimizing  solution  typically  yielded  excellent  model-data 
agreement:  scatter  plot  comparisons  (not  shown)  and  statistics 
(included  in  Table  4)  reflect  this. 

The  inferred  mud  thickness  values  for  the  three  zones  are 
shown  in  Fig.  9.  These  values  are  consistently  smaller  than  the 
observation-based  estimate  (40cm).  This  is  not  unexpected,  since 
in  the  blindfold  applications,  the  dissipation  is  overpredicted 
using  <5m0  =  40cm,  e.g.  a  -32cm  waveheight  bias  at  NDP  in  the 
one-dimensional  WDGL  set.  The  smaller  mud  thickness  given  by 
the  inversion  might  be  interpreted  to  mean  that  only  the  upper 
portion  of  the  mud  is  fluidized. 

The  Zone  B  solution  via  NDP  is  fairly  stationary  for  the  entire 
35-day  period,  at  <5mO  =  7-10crn,  while  the  Zone  C  and  D 
solutions  are  less  so;  for  example,  the  solutions  based  on  P5  are 
stationary  only  for  5-10  days  at  a  time.  There  is  similar  variability 
in  the  solutions  based  on  waveheights  observed  at  PA.  If  this 
variability  is  not  the  result  of  model  or  instrumentation  error,  it 
suggests  that  the  mud  layer  is  more  dynamic  than  the  infre¬ 
quently  sampled  data.  Such  variability  could  be  associated  with 
sediment  transport  or  time-variation  of  the  depth  of  fluidization 
of  the  mud  layer.  Further,  the  inversion  suggests  a  fairly  consistent 
non-zero  Sb0{  in  the  region  shoreward  of  the  9  m  contour,  where 
field  observations  do  not  suggest  the  presence  of  mud. 

An  assumption  discussed  above — that  the  waveheight  at  an 
observation  point  is  affected  only  by  the  mud  offshore  of  that 
point— appears  to  be  vindicated,  since  the  resulting  mud  solutions 
match  observed  waveheights  very  well.  For  example,  the  inversion 
solves  for  the  mud  thickness  in  Zone  B,  gets  a  match  to  observed 
waveheights  at  NDP.  progresses  to  and  creates  a  solution  for  Zone 
C.  and  this  modification  of  the  mud  thickness  in  Zone  C  does  not 
result  in  degraded  agreement  at  the  gage  at  the  deeper  end  of  this 
zone,  NDP. 


6.  Discussion 

6.  J.  Modeling  errors 

As  noted  already,  many  of  our  conclusions  require  that,  in  the 
simulations  with  Sbot  =  0,  this  omission  is  the  dominant  source  of 
error.  Further,  it  is  obvious  that  the  mud  thickness  values 
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Fig.  9.  Inferred  time  series  of  mud  thickness  in  cm.  This  shows  the  solution  in 
terms  of  with  <fWm>  held  constant. 
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determined  by  the  inversion  are  affected  by  all  errors  inherent  in 
the  modeling  and  measurements.  Thus,  it  is  useful  to  character¬ 
ize— and  if  possible,  estimate— this  error.  Model  error  is  caused  by 
error  in  forcing,  omissions  in  forcing,  errors  in  approximations  for 
physical  processes,  omission  of  physical  processes,  model  simpli¬ 
fications,  and  numerical  error.  Error  associated  with  uncertainty 
in  a  type  of  forcing — the  specification  of  the  thickness  and 
distribution  of  the  liquefied  mud  layer— was  addressed  in  Section  5. 
In  this  section,  other  types  of  errors  are  discussed. 

Omission  of  the  longshore  coordinate  in  the  inversion  is  one 
example  of  a  model  simplification  which  certainly  produces  error. 
However,  this  error  is  readily  estimated  from  comparison  of  Tables  3 
and  4  and  appears  to  be  modest.  Measurement  error  is  another 
concern;  however,  since  we  are  dealing  primarily  with  a  low- 
order  moment,  waveheight,  instrument  error  is  expected  to  be 
small:  less  than  15%  random  error  and  without  systematic  bias  as 
long  as  the  instrument  is  functioning  properly.  Error  in  bathy¬ 
metry  is  another  possible  source  of  error.  The  bathymetry  in  the 
nearshore  (depths  less  than  2  m)  can  be  expected  to  vary  in  time 
due  to  sediment  transport.  Jet-ski  surveys  (Holland  et  al„  this 
issue)  showed  no  appreciable  change  during  May  2005,  but  were 
unfortunately  not  conducted  duringjune  2005;  this  uncertainty  is 
relevant  to  model-data  comparisons  at  the  "P5”  location  during 
this  period. 

Spatially  varying  currents  cause  waves  to  refract  and  shoal/ 
strain  in  a  manner  similar  to  spatially  varying  water  depth. 
Further,  opposing  currents  can  cause  waves  to  steepen,  making 
them  more  likely  to  break.  Wave-current  interaction  is  not 
considered  in  the  presented  simulations.  Due  to  the  Patos  Lagoon 
inlet  only  5-10  km  to  the  northeast  of  Cassino  Beach,  surface 
currents  in  this  area  are  significant.  On  a  day  with  high  river 
discharge,  the  surface  currents  in  the  jet  from  this  inlet  can  easily 
exceed  1  m/s  (Vinzon  et  al„  this  volume).  If  this  jet  is  straight,  and 
assuming  linear  waves,  the  waves  will  be  refracted  inside  the  jet, 
but  will  leave  the  jet  with  their  original  characteristics,  and  thus 
currents  inside  the  jet  would  not  necessarily  affect  waveheights  at 
the  instrument  locations.  However,  the  currents  outside  the  jet 
are  also  quite  large,  sometimes  reaching  0.5  m/s,  with  a  recircula¬ 
tion  pattern  toward  the  northeast  observed  at  times.  A  typical 
group  velocity  of  a  wave  at  the  spectral  peak,  say  7=  10  s,  using 
15  m  water  depth,  is  9  m/s.  This  circulation  would  have  some 
effect  on  the  local  energy  level  of  such  a  wave  via  refraction  and 
straining,  but  not  a  large  effect.  For  shorter  waves,  the  group 
velocity  is  smaller,  e.g.  2.5  m/s,  and  the  effect  of  currents  could  be 
dominant;  however  since  these  waves  are  typically  in  the  spectral 
tail,  this  does  not  necessarily  imply  a  dominant  effect  on  the  total 
energy  (waveheight).  Therefore,  we  believe  that  currents  are  not  a 
major  source  of  error. 

Another  potential  source  of  error  is  the  assumption  of  uniform 
wave  conditions  along  the  offshore  boundary,  prescribed  equal  to 
the  spectra  observed  by  the  waverider  buoy.  Spatial  variability 
of  swell  fields  is  probably  not  a  problem,  since  swell  fields  tend 
to  be  much  larger  than  our  outer  grid  However,  since  our  off¬ 
shore  boundary  is  not  in  deep  water,  there  will  certainly  be 
some  variation  of  wave  conditions  along  the  boundary  associated 
with  shoaling,  refraction  and  non-conservative  processes  such 
as  bottom  friction  and  breaking.  Sensitivity  to  this  error  was 
estimated  by  running  a  regional  scale  model,  initialized  in  deep 
water,  with  output  along  the  offshore  boundary  of  our  sub¬ 
regional  scale  model.  This  test  was  performed  separately  for 
swells  approaching  from  the  south  (shore-oblique),  southeast 
(directly  onshore),  and  east-northeast  (shore-oblique).  In  the  case 
of  waves  approaching  from  the  south,  our  assumption  of  uniform 
boundary  conditions  is  estimated  to  result  in  20%  overprediction 
of  waveheight.  For  swells  from  the  southeast,  there  is  no  error 
(since  waves  at  the  nearshore  gages  are  coming  from  the  direction 


of  the  buoy),  and  for  waves  from  the  east-northeast,  a  small 
underprediction  (5%).  Comparing  bias  of  the  0  model  at  the 
NDP  location  versus  the  dominant  wave  direction,  there  is  not  a 
clear  correlation,  suggesting  that  this  effect  is  not  particularly 
important. 

Elastic  behavior  of  the  mud,  not  included  in  the  viscous  models 
used  here,  can  either  increase  dissipation,  especially  via  resonance 
(see  Zhang  and  Ng.  2006),  or  decrease  dissipation  by  restraining 
the  motion  of  the  bed  (see  MacPherson.  1980;  Hsiao  and  Shemdin, 
1980).  One  can  speculate  that  the  apparent  overprediction  of 
dissipation  by  either  Smiid  formulation  used  here  is  due  to  this 
omission.  Unfortunately,  without  actually  estimating  the  elasticity 
of  the  mud  at  Cassino  Beach  and  performing  calculations,  it  is  not 
possible  to  say  whether  including  this  process  would  increase  or 
decrease  estimated  dissipation.  Implementing  the  formula  of 
Zhang  and  Ng  (2006)  in  SWAN  is  a  potential  avenue  for  further 
work.  Treatment  of  the  mud  as  a  plastic  is  yet  another  option  (e.g. 
Mei  and  Liu,  1987).  In  any  event,  viscous  fluid  mud,  visco-elastic 
mud,  plastic  mud,  and  rigid  mud  may  have  coexisted  at  the  field 
site,  in  layers.  In  the  forward  model  applications,  it  is  assumed 
that  the  entire  mud  layer  is  liquefied.  As  shown  in  Section  3.2, 
if  only  a  small  portion  of  the  mud  layer,  say  the  upper  10  cm, 
is  liquefied,  this  would  tend  to  reduce  dissipation.  Thus,  this  is 
another  possible  explanation  for  overpredicted  dissipation  with 
the  forward  models. 

Due  to  lack  of  knowledge  regarding  sand  bedforms  on  this 
beach,  bottom  friction  is  not  applied  in  the  simulation  in  areas 
where  mud  thickness  is  zero.  This  omission  could  account  for 
some  of  the  underprediction  of  dissipation  shoreward  of  the  9  m 
contour  suggested  by  the  inverse  modeling,  but  we  do  not  believe 
that  the  bedforms— if  they  existed  at  all  were  large  enough  to 
play  a  major  role.  This  is  discussed  quantitatively  in  Section  6.3. 

6.2.  Relative  importance  of  depth  limited  breaking 

Though  not  included  here,  comparisons  were  made  to  a  SWAN 
model  implementation  with  depth-limited  breaking  disabled. 
These  model  results  suggest  surf  breaking  is  only  significant 
shoreward  of  the  PA  location  and  only  during  the  higher  wave 
events.  This  is  consistent  with  traditional  rules  of  thumb 
regarding  maximum  height  to  depth  ratios  (e.g.  Dean  and 
Dalrymple,  1991).  Thus,  breaking  is  only  relevant  to  the  compar¬ 
isons  at  the  P5  location,  and  then  only  for  a  few  time  periods.  The 
comparison  suggests  that  during  the  June  16-17  high  wave  event, 
surf  breaking  plays  a  dominant  role  in  the  P5  comparisons.  In  fact, 
at  this  location  the  highest  waves  do  not  occur  during  June  16-17, 
but  during  periods  prior  and  after,  June  14  and  June  20  when 
breaking  was  apparently  reduced  Encouragingly,  this  behavior  is 
seen  both  in  the  observations  and  in  the  model  that  includes  surf 
breaking. 

6.3.  Application  with  sandy  bed 

It  was  determined  in  Section  4  that  inclusion  of  dissipation  by 
wave-mud  interactions  S ^  greatly  improved  bias  statistics, 
especially  at  gages  PA  and  P5.  The  question  might  be  asked, 
"Could  the  same  improvement  be  achieved  by  using  a  formula  for 
Sbot  based  on  a  sandy  seabed?"  Such  a  study  was  in  fact 
performed.  Specifically  171  two-dimensional  simulations  from 
Section  4  were  applied  using  the  Madsen  et  al.  (1988)  formula  for 
dissipation  by  bedforms.  Since  there  was  no  observation-based 
guidance,  a  bedform  amplitude  of  5  cm  was  applied  uniformly 
over  the  domain.  Interestingly,  these  results  show  a  similar 
improvement  in  bias  compared  to  the  case  with  Srnud.  Comparing 
the  statistics  with  mud  dissipation  versus  those  with  sand 
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dissipation,  they  were  fairly  similar,  with  the  former  being 
moderately  better  at  P5  and  the  latter  slightly  better  at  NDP  and 
PA.  However,  the  similarity  of  magnitude  of  the  mud-based  and 
sand-based  dissipation  is  only  due  to  the  application  of  the  sandy 
bedforms  over  a  larger  area. 

The  actual  decay  rate  based  on  Madsen  et  al.  (1988)  with  5  cm 
bedform  amplitude  is  considerably  smaller  than  that  of  Ng  as 
applied  in  Section  4.  Unfortunately,  comparison  of  the  decay  rate 
based  on  Madsen  et  al.  (1988)  with  those  in  Table  1  is  not 
straightforward,  because  this  source  term  is  nonlinear,  meaning 
that  the  exponential  decay  rate  of  the  wave  amplitude  is 
dependent  on  the  amplitude  itself.  But  to  provide  an  example, 
in  the  case  of  h  =  12  m  in  Table  1,  if  we  use  a  waveheight  of  1  m 
and  a  bedform  amplitude  of  5  cm,  the  decay  rate  k,  =  3.5  x  10  5 
rad/m  compared  to  k,  =  1.3x10  4  rad/m  using  Ng  with  stated 
mud  characteristics. 

6.4  /nversion  for  oddirionol  voriables 

This  inversion  described  above  assumes  constant,  uniform 
values  for  the  mud  density  and  viscosity,  pm  =  1310  kg/m3,  and 
vm  =  7.6x10  3m2/s,  both  based  on  rheometry.  Only  the  mud 
thickness  and— in  a  limited  fashion— the  horizontal  (x)  distribu¬ 
tion  are  unknowns  to  be  solved  for.  Thus,  the  mud  thickness  being 
solved  is  <5m<0  =  ^m.o(xT)  and  the  outcome  of  the  inversion  is  the 
solution  set  in  terms  of  <t5m.oA*t>-  However,  as  noted  in  Section 
3.2,  there  is  also  a  range  of  uncertainty  in  the  mud  viscosity  and 
density,  so  the  solution  set  might  be  more  generally  in  terms 
of  <  > •  An  obvious  question  is,  "If  we  had  used  a 

different  mud  viscosity,  within  the  range  of  probable  values, 
would  we  have  recovered  mud  thickness  values  closer  to  the 
observation-based  estimate  (40cm),  instead  of  20cm?”.  To  repeat 
the  inversion  process  for  multiple  possible  values  of  pm  and 
vm — say  ten  values  each  would  increase  the  total  computational 
requirement  by  a  factor  of  100.  However,  this  is  not  strictly 
necessary.  Instead,  for  each  <x,f>,  we  can  use  the  WDGL 
formulation  to  calculate  k,  from  <$m0  given  the  <vm,pm>  values 
used  in  the  actual  inversion  and  then  apply  the  formulation  again 
to  calculate  combinations  of  that  produce  this  k,% 

None  of  this  requires  new  computations  with  the  wave  model. 
Application  of  the  WDGL  formulation  in  this  manner  does  require 
a  representative  value  of  wave  period  T,  so  it  is  only  an 
approximate  solution,  but  the  outcome  appears  to  be  only  weakly 
sensitive  to  variations  in  the  chosen  values  for  T. 

7.  Summary  and  conclusion 

Dissipation  by  viscous  mud  by  two  formulations  was  inde¬ 
pendently  implemented  in  SWAN:  (1)  De  Wit  (1995)  in  Winter- 
werp  et  al.  (2007)  and  (2)  Ng  (2000)  herein.  Both  are  applied  to 
the  Cassino  Beach  2005  field  experiment.  Calculations  with  the  Ng 
formula  reveal  that — though  this  is  a  clearly  muddy  area 
dissipation  by  mud  in  this  field  experiment  is  weak  relative  to 
examples  given  in  the  literature,  as  well  as  a  hypothetical  scenario 
where  the  mud  is  at/near  the  surf  zone  in  Cassino,  such  as  occurs 
there  periodically.  Calculations  also  suggest  that  the  uncertainty 
in  the  specification  of  the  mud  in  the  modeling  exercise  is  likely  to 
have  significant  impact  on  results.  Uncertainty  in  viscosity  and 
fluidized  mud  thickness  are  found  to  be  particularly  important. 
Further,  the  strong  sensitivity  of  dissipation  to  the  local  water 
depth  implies  that  the  modeling  is  highly  sensitive  to  errors  in  the 
horizontal  distribution  of  the  mud. 

The  two  formulations  are  compared  herein.  There  are  two 
primary  differences.  First,  the  Ng  formula  is  intended  only  for 
cases  where  the  mud  layer  is  thin  relative  to  the  overlying  water 


layer,  and  it  is  shown  here  that  the  WDGL  implementation  has  a 
broader  range  of  validity.  However,  the  thickness  of  the  mud  in 
the  Cassino  Beach  application  does  appear  to  be  well  within  the 
limits  of  validity  of  the  Ng  formula.  The  second  major  difference  is 
the  method  of  converting  the  exponential  decay  rate  k,  into  a 
spectral  dissipation  term  Sbot:  the  methods  predict  similar  k,  for 
the  case  study,  but  Sbot  is  significantly  larger  with  the  Winterwerp 
et  al.  (2007)  method. 

In  the  Cassino  Beach  2005  application  of  the  wave  model  with 
Sbot  =  0  (i.e.  no  dissipation  by  wave-bottom  interaction),  wave 
energy  is  overpredicted  at  all  three  measurement  locations. 
Assuming  that  other  sources  of  systematic  error  are  small,  this 
overprediction  of  energy  suggests  an  underprediction  of  dissipa¬ 
tion  This  trend  is  quite  consistent  for  the  35-day  period  studied, 
with  the  only  exception  being  during  a  single  high  energy  wave 
event.  During  this  time  period,  the  waveheight  is  well  predicted 
using  3^  =  0.  Two  possible  explanations  are  suggested:  (1)  the 
water/mud  interface  may  have  been  obliterated  by  sediment 
suspension  during  the  high  wave  event,  or  (2)  the  mud  may  be 
thixotropic,  i.e  reduced  viscosity  under  greater  wave  forcing. 

At  all  other  time  periods,  application  of  S^t  =  Smud,  with  mud 
thickness,  horizontal  distribution  and  rheology  based  on  observa¬ 
tions,  results  in  modestly  underpredicted  wave  energy  for  both 
dissipation  formulations  applied,  suggesting  that  the  dissipation 
may  be  overpredicted.  We  offer  three  speculative  explanations  for 
this.  First,  the  utilized  models  assume  zero  elasticity;  we  point  out 
that  the  use  of  non  zero  elasticity  might  improve  results  by 
reducing  the  predicted  dissipation  (e.g.  see  MacPherson,  1980), 
Second,  the  actual  mud  was  certainly  less  uniform  than  as  applied 
in  the  model;  thus,  it  might  have  been  significantly  thinner  in 
places.  Third,  there  is  a  very  strong  possibility  that  only  the  top 
portion  of  measured  mud  was  fluidized. 

The  forward  modeling  methodology  described  above  is 
potentially  useful  in  wave  hindcasting  and  forecasting,  if  the 
rheology  and  mud  distribution  are  known  with  some  degree 
of  confidence.  In  the  reverse  situation,  the  waves  are  measured 
and  the  mud  is  poorly  described.  A  method  of  inferring  mud 
distribution  from  observed  wave  height  distributions  is  presented. 
In  application  of  this  method,  the  observed  wave  heights  are 
successfully  recovered  by  the  inversion.  The  optimal  thicknesses 
showed  both  spatial  and  temporal  variability  suggesting  that 
more  detailed  observations  of  these  difficult  parameters  are 
necessary  to  properly  validate  dissipation  mechanisms. 
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Appendix  A.  Wind  vector  algorithm 

Two  sources  were  initially  available  to  provide  time  series 
of  wind  vectors  for  the  time  period  of  the  field  experiment: 
an  anemometer  located  near  Cassino  Beach,  and  the  Navy’s 
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operational  global  atmospheric  model,  NOGAPS  (Hogan  and 
Rosmond,  1991).  Comparing  these  two,  the  wind  directions  were 
consistent,  but  the  anemometer  wind  speeds  and  the  NOGAPS 
winds  speeds  over  a  nearby  land  point  (53  W,  32°S)  were 
systematically  lower  than  the  NOGAPS  wind  speeds  over  a  nearby 
sea  point  (52  W.  33  S).  presumably  due  to  land  frictional  effects. 
Direct  use  of  the  anemometer  winds  was  ruled  out  due  to  this 
apparent  bias.  There  was  concern  that  temporal  structure  of  the 
NOGAPS  wind  speeds —available  only  every  three  hours— do  not 
correspond  with  that  of  the  boundary  forcing,  and  further  may 
contain  phase  errors.  Thus,  a  third  independent  estimate  of  wind 
vector  was  derived  based  on  observed  high-frequency  wave 
energy.  An  automated  procedure  works  as  follows. 

•  A  sequence  of  SWAN  simulations  were  performed,  each  with  a 
different  10  m  wind  speed,  ranging  from  4  to  25  m/s,  at  1  m/s 
increments.  The  grid  design  was  similar  to  that  used  in  the 
actual  hindcasts,  but  there  was  no  boundary  forcing  and  the 
wind  direction  was  always  from  318  (directed  offshore)  These 
simulations  were  used  to  create  a  database  on  variation  of  one- 
dimensional  spectra  with  wind  speed. 

•  For  each  wavender  spectrum,  the  high-frequency  tail  is 
identified. 

•  The  energy  level  in  this  frequency  range  (specifically,  the 
frequency- integrated  variance)  is  compared  with  a  similar 
integration  of  the  database  spectra  and  the  most  closely 
matching  wind  speed  is  selected  for  use  in  the  actual  hindcast. 

•  The  wind  direction  is  assumed  to  be  identical  to  the  mean 
direction  of  the  high-frequency  tail. 

The  measured  spectra  were  often  of  complex  mixed  sea/swell 
conditions,  so  the  most  significant  challenge  with  this  automated 
procedure  was  to  correctly  identify  the  windsea  portion  of  each 
spectrum  without  a  priori  knowledge  of  the  wind  speed.  A 
quality-control  procedure  was  developed  to  flag  dubious  solutions 
and  exclude  these  cases  from  the  set  of  hindcasts. 

The  inferred  winds  are  specific  to  a  single  point,  the  wavender 
buoy  location.  In  the  wave  model  applications,  this  is  applied 
uniformly  over  the  domain,  so  the  reduction  of  wind  speeds  by 
frictional  effects  over  land  is  not  accounted  for  in  the  wave  model. 

The  three  independent  estimates  of  wind  speed  are  shown  in 
Fig  Al;  here  the  NOGAPS  time  series  is  for  the  sea  point  (52  W. 
33  S).  The  mean  wind  speed  of  the  inferred  and  NOGAPS 
estimates  is  roughly  similar,  and  the  anemometer  estimates  are 
much  smaller  than  both.  The  inferred  and  NOGAPS  wind  speeds 
are  expected  to  be  more  representative  of  the  winds  that  would 
generate  waves  arriving  at  the  nearshore  gages.  Selection  of  one  of 
the  two  is  subjective,  since  each  contained  unknown  errors. 
NOGAPS,  since  it  is  a  global  model,  uses  relatively  coarse  native 
geographic  resolution,  which  is  potentially  a  problem,  especially 
near  coastlines.  The  inferred  wind  speeds  also  contain  errors, 
being  sensitive  to  errors  in  the  wave  generation  physics  of  the 
wave  model.  However,  since  these  errors  in  wave  generation 
physics  would  be  replicated  in  the  actual  hindcasts,  it  can  be 
argued  that — in  the  hindcast  application— the  inferred  wind 
speeds  are  more  likely  to  produce  windsea  growth  similar  to  that 
which  is  observed,  even  though  the  wind  speeds  themselves  are 
not  necessarily  more  accurate  than  those  of  NOGAPS.  Further,  the 
temporal  structure  of  the  inferred  wind  speeds  corresponds 
directly  with  that  of  the  boundary  forcing,  in  contrast  to  NOGAPS. 
Thus,  the  inferred  wind  speeds  are  selected  for  the  hindcasts. 

As  was  mentioned  in  Section  4,  97  of  528  time  periods  were 
discarded  due  to  mismatch  of  wave  height  at  the  offshore 
boundary,  indicative  of  problems  with  wave  growth  internally. 
This  was  a  surprising  result,  since  the  wind  speed  inference 


Fig.  A1.  Time  senes  of  wind  speed  by  three  independent  estimates. 


algorithm  should  provide  the  optimal  wind  speed  for  producing 
in  the  model  the  same  wind  sea  as  observed  at  the  buoy. 
The  discrepancy  was  found  to  be  due  to  non-physical  dependence 
of  the  whitecapping  term  on  mean  steepness:  a  problem  that  is 
improved — but  apparently  not  eliminated— using  the  van  der 
Westhuysen  et  al.  (2007)  source  terms.  Specifically,  windsea 
component  propagating  nearly  parallel  with  the  offshore  bound¬ 
ary  was  found  to  grow  too  quickly  when  swell  is  present.  Since  the 
wind  speed  inference  algorithm  is  applied  without  swell  present, 
inferred  wind  speeds  are  not  affected  by  this  problem,  while  the 
actual  simulations  are.  In  any  event,  the  decision  to  omit  these 
simulations  was  somewhat  subjective,  since  it  was  not  actually 
determined  whether  this  spurious  windsea  component  signifi¬ 
cantly  affects  nearshore  comparisons. 
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